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Abstract 

In this paper, we consider the problem of multi-resolution compressed sensing (MR-CS) reconstruc¬ 
tion, which has received little attention in the literature. Instead of always reconstructing the signal at 
the original high resolution (HR), we enable the reconstruction of a low-resolution (LR) signal when 
there are not enough CS samples to recover a HR signal. We propose an approximate message passing 
(AMP)-based framework dubbed MR-AMP, and derive its state evolution, phase transition, and noise 
sensitivity, which show that in addition to reduced complexity, our method can recover a LR signal 
with bounded noise sensitivity even when the noise sensitivity of the conventional HR reconstruction 
is unbounded. We then apply the MR-AMP to image reconstruction using either soft-thresholding or 
total variation denoiser, and develop three pairs of up-/down-sampling operators in transform or spatial 
domain. The performance of the proposed scheme is demonstrated by both ID synthetic data and 2D 
images. 
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I. Introduction 


Recently compressed sensing (CS) has been studied extensively as an efficient way of acquiring and 
reconstructing sparse signals 0 - Many CS reconstruction algorithms have been developed, e.g., convex 
optimization Q, greedy method Q, iterative thresholding 0, and approximate message passing (AMP) 

0 -®. 


The AMP is a particularly attractive framework, due to its near-optimal reconstruction performance, 
low complexity, and the capability of predicting its performance from its state evolution. This leads to 
the discovery of the phase transition property of the AMP, which states that when the sampling rate of a 
sparse signal is below a threshold defined by a phase transition curve (PTC) 0, 0, | |T0[ , | [TT[ , the CS 
algorithm will fail to recover the signal with high probability even if there is no sampling noise. In the 
noisy case, the noise sensitivity is unbounded, where the noise sensitivity is the minimax mean squared 
error (MSB) of the reconstruction. This is analogous to the rate-distortion bound in information theory. 
Therefore, in applications that a large amount of CS samples need to be transmitted to a receiver, the 
receiver has to wait until it receives enough samples before it can recover the signal. This could incur 
undesired delays. 

This paper is motivated by the following fundamental question: if in the case above we are allowed 
to reconstruct low-resolution (LR) previews instead of the original high resolution (HR) signal, can we 
recover high-quality LR signals so that we can enlarge the feasible operating region of the system? We 
call this framework CS with multi-resolution reconstructions, or MR-CS for short. It opens up many 
questions. For example, how to design the sampling and reconstruction algorithms? What is the highest 
resolution that can be reconstructed at each sampling rate? What are the expressions of the phase transition 
curves for different LR reconstructions? In addition, a straightforward approach is to first reconstruct a 
HR signal using existing reconstruction methods, and then downsample it. Therefore another question is 
how much gain we can get over this simple method. Note that a carefully designed LR reconstruction 
algorithm should at least have lower complexity than this simple method, because it could reconstruct 
the LR signal directly. 

Although the need for multi-resolution (MR) or scalable reconstruction has been well recognized in 
multimedia transmission, leading to the development of standards such as JPBG 2000 and H.264/SVC 
|T2), I®, the problem has received little attention in CS. The schemes that are most relevant to ours 
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are In p4} , some rules are proposed to design efficient up-/down-sampling matrices for MR 

reconstruction, and the number of nonzero entries of the LR image in transform domain is shown to be no 
larger than that of the HR image. Therefore the required sampling rate for stable LR image reconstruction 


is less than that of the HR reconstruction. However, the analysis in 1141 is qualitative, and only some 
loose bounds are provided. Moreover, the impact of the MR design on the quality of the measurement 
matrix is not studied, which can be measured by, e.g., restricted isometry property (RIP) constant Q 
and mutual coherence Q. Besides, only the noiseless case is considered in p4|. 


A similar problem to ours is studied in |16|, where two solutions are proposed. In the first method, 
the sampling matrix is designed to have non-uniform sampling, which is quite restrictive, since the 
matrix should be redesigned whenever a new result with different resolution is needed. The second 
method modifies the sampled data of the HR image to be close to the data acquired directly from the 
target LR image. Although it works empirically, there is no theoretical guarantee. In addition, although 
it is mentioned in m that the CS sampling rate for the LR reconstruction is increased, the change 
of the sparsity rate is not considered. Moreover, the complexity of this approach is even higher than 


reconstructing the HR image directly. We will show in this paper that the second solution in |16| is a 
special case of our proposed MR-AMP framework. 


Recently, a special two-resolution CS reconstruction scheme is proposed in |15|, where the sampling 
matrix is designed such that a LR reconstruction can be obtained by direct matrix inversion. 

The MR concept has also been used in some CS schemes such as with different purposes 

from ours. In p7|, Bayesian CS is used to detect the primary user in cognitive radio. It first performs 


the detection in LR, and then refines the signal around the detected primary user spectrum. In 1I8|, a 
CS-based two-layer scalable image coding is proposed, where the encoder employs two measurement 


matrices with different sizes, and inter-layer prediction is used to reduce the bit rate. In 1191, the authors 


extended the Kronecker CS |22| to MR measurements, such that the sensing is performed on the LR 


image, and the goal is to recover the HR signal from LR measurements. In |20|, a multiscale framework 
is proposed for CS of videos. The motion vectors are estimated at different resolutions and served as 
the input to higher resolution frame recovery. The sensing is applied to different resolutions of the same 
frame. In our proposed framework, the sensing is only performed on the original HR image. Therefore, 


the framework in |20| is more like source coding, but not sensing and coding simultaneously. In |2I |, the 
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authors use advanced denoising methods in the multiscale wavelet domain to improve the performance 


of AMP reconstruction, similar to | [T()| . However, the reconstruction still has the same resolution as the 
source. 

In this paper, we develop a general theory for MR-CS reconstruction, and propose a MR-AMP algorithm 
to reconstruct a LR signal if the sampling rate is too low. Our method does not impose any constraint on 
the measurement matrix. Therefore it enables more LR reconstruction choices. Also, theoretical analysis 
can still he obtained. Instead of having only one phase transition curve (PTC), we obtain a family of 
PTCs that specify the sampling rate thresholds to get bounded noise sensitivity with different resolutions. 
Moreover, the noise sensitivity is derived explicitly. The performance of the proposed scheme is verified 
using both synthetic data and natural images. 

The rest of this paper is structured as follows: Sec. |T^ presents the mathematical model of MR-CS 
problem and provides the necessary conditions that the MR up/down-sampling matrices should satisfy. 


Sec. Ill is devoted to MR-AMP algorithm and its updating rule. Sec. IV establishes the theoretical 
analysis of MR-AMP. Sec. |V] discusses the application of MR-AMP in images, and develops three sets 
of up/down-sampling matrices. Sec. |V^ presents simulation results, validates the state evolution of MR- 
AMP, and gives guidelines on tuning the parameters of the algorithm. It also compares the performance 
of MR-AMP with the original HR-AMP with different denoisers, in terms of reconstruction quality and 


algorithm complexity. Some preliminary results of this paper are reported in |23|. 


II. Formulation and Conditions of MR-CS Reconstruction 
The goal of the classical CS is to recover a ni x 1 vector x from a m x 1 noisy measurement y with 
m < ni, i.e., 

y = Ax + w. (1) 

In this paper, entries of the m x ni measurement matrix A are i.i.d. Gaussian with zero mean and a 
variance of 1/m, denoted by AA(0,1/m). Each entry of the noise vector w also follows i.i.d. Gaussian 
distribution with zero mean and a variance of The CS undersampling ratio is defined as (5i = m/ni. 

Since the system is underdetermined, it cannot be solved without exploiting the special structure of x. 
Some examples of structured signals are given in Q, including simple sparse signals, block sparse signals, 
mostly constant non-decreasing signals, and piecewise constant signals. Following the notations in Q, 
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the family of probability distributions for a particular type of structured signals over is denoted as 
where ei < 1 is a constant sparsity ratio, and the expected amount of useful structured information 
in the signals is at most ki = niCi- The definition of the useful structured information depends on the 
nature of the structure. Let denote a distribution in and x be a signal with distribution Vm. 

In this paper, we focus on the following two families of structured sparsity. 

Definition 2.1: The family of distributions that generates simple sparse signals is defined as (Eq. (1.2) 
in @) 

{||x||o} ^ niEi] , (2) 

where fhe Iq norm ||x||q denotes the number of nonzero entries of vector x. Therefore, the expected 
number of non-zero entries of signals in this family is at most niSi. 

Definition 2.2: The family distributions that generates piecewise constant signals is defined as (Sec. 
V in |§) 

-pPC ^ 

ni,£i — 

(3) 

{uni : {#{f G [l,ni - 1] : Xt+i Xt}} ^ uiEi} , 

where #{•} denotes the number of times the condition in the operator is true. Therefore, the expected 
number of change points within signals of this family is at most niEi. 

In the proposed MR-CS reconstruction framework, instead of always recovering the signal with the 
original resolution ni, we allow the reconstruction of various lower resolution signals Ud (nd < ni) when 
the number of available CS samples is too small. 

The MR downsampling factor is defined as 


d = ni/nd. (4) 

Note that this MR downsampling factor should not be confused with the CS undersampling ratio 
(5i = m(ni. In this paper, we are interested in the case m < Ud, i.e., the recovery of the LR signal is 
still an underdetermined CS problem. The equivalent CS undersampling ratio for the LR reconstruction 
is 6d = rn/ud = d6i > (ii. Let kd be the expected amount of useful information contained in the 
LR signal. The expected sparsity ratio of the LR signal is Ed = kd/ud. We also define another factor 
Pd = ^dldd = kd/m. Clearly, a signal with larger pd needs more measurements (larger 5d) to recover. 
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Let Drf be a X ni downsampling matrix, be a ni x upsampling matrix, and = D^x be 


the Ud X 1 downsampled version of x. The LR-CS problem can be formulated as 114| 


y = Ax + w = A(UrfXrf + X - UdXrf) + w 
= AUrfXrf + A(I - UdDrf)x + w, 


(5) 


where AU^ is the equivalent measurement matrix for the LR signal x^^, and A(I — UdDrf)x is the 
additional approximation error term when x^ is the target signal to be recovered. Note that this error 
term depends on the signal x. 

The downsampling and upsampling matrices and play an important role in the MR-CS. In this 
paper, we require them to satisfy three conditions. 

Condition 2.1: The downsampling and upsampling matrices and should be chosen such that 
if we first upsample a LR signal and then downsample to the original resolution, we can get back the 
original LR signal without any error. That is. 




( 6 ) 


Since is a tall matrix, this mild condition can be easily satisfied. In 115|, fhe authors design a 


special two-resolution CS system such that a m x 1 LR signal can be recovered directly from the m x 1 
CS sample y. This can be considered as a special case of our setup. 

The second condition is about the quality of the measurement matrix for the LR reconstruction. 

Condition 2.2: The quality of the equivalent measurement matrix for the LR reconstruction should be 
no worse than that of the HR reconstruction. 

For different reconstruction algorithms, different criterions are used to evaluate the quality of the 
measurement matrix, e.g., the RIP constant for the basis pursuit algorithm Q and the mutual coherence 
for the orthogonal matching pursuit algorithm Q. The solution in this paper is based on the AMP 
algorithm; hence we follow the requirement in Q, Q, p0| , pT| that each entry of the LR measurement 
matrix should be i.i.d. Gaussian with zero mean and a variance of 1/m. 


Since m < in our case, the MR-CS problem here cannot be solved directly without exploiting the 
structure of x^. Moreover, the LR signal should be easier to recover than the HR signal, i.e., the amount 
of useful information kd contained in x^ should be no more than the amount ki in the original HR signal 
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X. We therefore also require the downsampling matrix to satisfy the following condition. 

Condition 2.3: If x belongs to the family iFni,ei in basis the downsampling matrix should be 
chosen such that = D^x belongs to the family iFnd,Sd in basis = {D^^} — {0} with Ed ^ dci- 


Some results similar to Cond. 2.3 have been reported in |14| for simple sparse vectors, which is a 


special case of Cond. 2.3 as summarized below. 

Condition 2.4: If x is sparse in basis 'J', then x^ = D^x is sparse in the non-zero projected low- 
dimension basis = {D^^} — {0}. The sparsity kd of x^^ is no larger than k, the sparsity of x, if the 
columns of are linearly independent. 


Our condition in Cond. 2.3 is not restricted to simple sparse vectors, and can be used for other special 
structures that x follows, such as piecewise constancy. 

In Sec. 1^ we will design three pairs of up-/down-sampling matrices for images that satisfy the three 
conditions above perfectly or approximately. One pair is for simple sparse signals and two pairs are 
for piecewise constant signals. The conditions listed above can also be used to design matrices for the 
multi-resolution reconstructions of other types of structured sparse signals. 

Note that the term ’’multi-resolution” in our paper is slightly different from that in the wavelet transform 
literature, because our method only reconstructs each of these LR signals independently, and how to use a 
LR reconstruction to help a HR reconstruction is not addressed in this paper. Nevertheless, we will show 
in Table |VII| that by simply upsampling the recovered LR image to the target HR, we can sometimes 
provide better HR image than reconstructing the HR image directly from the measurements. 


HI. Multi-Resolution Approximate Message Passing 

In this section, we propose an approximate message passing (AMP)-based algorithm to solve the MR- 
CS problem. Without loss of generality, we assume that the signal belongs to the structured sparse family 
Tni,ei in the canonical basis. 

The main idea of the original AMP is to transform the CS reconstruction problem into a denoising 
problem Q, i.e., estimating Xo from its noisy observations Xod-cre, where entries of e are i.i.d. Gaussian 
with zero mean and unit variance, and u is a constant. In each iteration of AMP, a pseudo-data z* = 
X* -I- is first formed. It is then denoised by a denoising function r), where a* is the standard 

deviation (std) of z* and r is the tuning parameter of the denoiser. Finally the residual of the measurements 
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X 


■t+i _ 



(7) 


pt+i = y - Ax*+^ + 6‘r* 


where is the Onsager term, which is related to the divergence of the denoiser hy 


6 ' 


m 


1 


div7/^t-i(u;r)| 



( 8 ) 


For different structured signals, different denoisers should he used. For example, for simple 

sparse signals, the well-known soft-thresholding should he used, whereas total variation (TV) denoiser 
is more appropriate for piecewise constant signals 1^. 

In order to apply AMP to the MR-CS problem in Eq. Q, we propose the following multi-resolution 
approximate message passing algorithm (MR-AMP), 



(9) 


where A^ = AU^A is the corresponding measurement matrix for the LR reconstruction, with A being a 


diagonal matrix determined by the upsampling matrix to normalize the colunms of AU^. 6^ is similar 


to Eq. (j^ except that m becomes Ud- Instead of estimating Xq, we are trying to estimate Xd,o = D^Xq 
from the pseudo-data with std cr^ using the denoising function 

The original AMP in Eq. (j7]| is a special case of MR-AMP in Eq. Q with d = 1. In this paper, 
we denote the original AMP as high-resolution approximate message passing (HR-AMP), and MR-AMP 
with d > 1 as low-resolution approximate message passing (ER-AMP). Since the dimensions of A^ and 
x^ are smaller than those of A and x, the complexity of ER-AMP is thus lower than HR-AMP. Note 
that the proposed MR-AMP does not impose any additional constraint to the measuring matrix A in the 
original AMP. It only modifies the reconstruction algorithm to get different ER estimates of the signal. 
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IV. State Evolution and Phase Transition of MR-AMP 

In this section, we analyze the theoretical performance of the proposed MR-AMP in terms of its state 
evolution, phase transition, and noise sensitivity. 


A. State Evolution 

The availability of the state evolution analysis is an important advantage of AMP over many other 
CS algorithms. Empirical findings show that the MSEs of AMP with various denoisers can he predicted 
accurately hy its state evolution Q, |10|, which descrihes the asymptotic limit of the AMP estimates 
in Eq. (j^ when m, ni —)> oo, for any fixed t | |TT| . Starting from 9^ = ||xo|| 2 /ni, the state evolution 
generates a sequence of numbers through the following iterations. 

1 

( 10 ) 


= ^6i*(xo,(5i,c7^,r) + f7^, 


1 


6»*+^(xo,(5i,(T^,r) = —E ||??a*(xo + cr*e;r) - XqI 


m 


where the expectation is with respect to e ~ A((0,I). Eor large values of m and ni, the state evolution 
predicts the MSE of the AMP algorithm in Eq. §• i.e., 0*(xo,(ii,(T2,r) ss ^ ||x* -XqH^. 

To get the state evolution of the proposed MR-AMP, we start from 9^ = ||X(^ where x^^ o is the 

target ER signal. Eet ci^ ^ denote the variance of the MR-AMP noise in Eq. (j^, including contributions 
from the approximation error and measurement noise, which is equal to -|- 1/m || (I — UrfD(i)x|| 2 ), as 


will be shown in Sec. IV-C The state evolution of the MR-AMP is thus given by the following iterations. 

^d, Tj -I- (72^^, 

( 11 ) 

12 ’ 


(<^df = ^^d(xd,o, ^d, r) + 


9*/\xd,o,^d,(7l^,T) = —E||r/^t(xrf,o + o-de;r) - Xd,of 

7td 


where cr^ is the predicted std of the estimate in Eq. d9|. If d = 1, Eq. (Ill reduces to that of AMP 


in Eq. ([1^). 

Note that the state evolution of AMP is only proved rigorously for scalar denoisers, but not for non¬ 


scalar denoisers, such as total-variation-based denoisers and other more advanced denoisers |10|, |21|, 


|. However, similar to observations in these papers, empirical findings in Sec. |Vl] show that in all 
cases studied in this paper the MSEs of the MR-AMP can be predicted accurately by the state evolution 


above. 
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B. Noiseless Phase Transition of LR-AMP 


In CS reconstruction without sampling noise, the phase transition curve (PTC) defines the minimum 
number of CS measurements required to perfectly recover Xq, i.e., 0°°(xo, 5i, 0, r) —)■ 0 |[^. In this 
part, we investigate the noiseless phase transition of MR-AMP, where we assume both cr^ = 0 and 
II A(I — UdDrf)x ||2 = 0 in Eq. The latter is possible for some special signals, and an example will 


be given in Sec. VI We will show that by allowing LR reconstruction, the MR-AMP admits a family 
of PTCs, thereby enabling perfect reconstruction of a LR signal in the infeasible region of the original 
HR-AMP This is an important generalization of the AMP theory. 

The family P'n,s is scale-invariant |[^, i.e., rjaiy, r) = cJ7?i(y/(T; r). Therefore we only need to consider 
a = 1, and we can simplify the notation r]a{y;T) as r/(y;r). We then define the following asymptotic 
minimax MSE when a denoiser rj with parameter r is used to recover signals in the structured sparse 
family Tnua @- 


1 9 

M{8i\r])= lim —inf sup ||r/(xo-f e; r) - XqH^ , (12) 

ni^ooni r 

In words, M{si\r]) is obtained by tuning the denoiser parameter to minimize the MSE per coordinate of 
the least favorable distribution in the family. The tuning rules of the parameters r are provided in Sec. 

NUB 

The minimax MSE has some basic properties Q, Q. Eirst, since the denoising can improve the 
reconstruction, we have 0 < M{si\rj) < 1. Besides, M{si\r]) —)• 0 when ei —)• 0, and M{ei\r]) —)• 1 when 
ei —)• 1. Second, M{si\r]) is monotonically increasing with respect to ei Q, because the reconstruction 
difficulty increases with si. 

The detailed expression of M{ei\ri) for AMP with various denoisers is derived in Q, |[^, pT| . More 
importantly, it is shown in Q that M{ei\r]) defines the minimum CS undersampling ratio (5i for perfect 
reconstruction, i.e., it describes the phase transition curve of AMP as follows. 

Theorem 4.1: In the noiseless case, when using AMP with denoiser t] to reconstruct signals in 
the AMP succeeds with high probability if 


> M{ei\r]). 


( 13 ) 


Viceversa AMP fails with high probability for Ji < M{£i\r]). 
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Combining Theorem 4.1 and the conditions in Sec. |T^ we obtain the following generalized phase 
transition result for MR-AMP, which specifies the minimum sampling ratio to perfectly recover a LR 


signal. When d = 1, it reduces to Theorem 4.1 


Corollary 4.2: When Cond. 2.1 2.2| and [23] are satisfied, if a signal x G is sampled according 


to Eq. and if = 0 and ||(I — UdD(i)x ||2 = 0 in Eq. ([^, then a ER signal x^ G with 

Ed ^ d El can be reconstructed perfectly with high probability via the ER-AMP in Eq. Q when the CS 
undersampling ratio satisfies 

5i> M{dEi\rf)/d, (14) 

where M(ei|r?) is the minimax MSE of the original HR-AMP. Viceversa the ER-AMP fails with high 
probability for di < M{dEi\'q)/d. 

Proof: As mentioned before, 6d = d5i. Since there is no approximation error in Eq. (|^, Theorem 


4.1 can be applied directly to the ER-AMP. Therefore the ER-AMP succeeds with high probability if the 
CS sampling ratio satisfies 

5d = d5i > M{Ed\ri). 

If Cond. |2.3| is satisfied, we have Ed ^ dci- Eq. ( [14] ) can thus be obtained using the property that 
M{Ed\il) is monotonically increasing with respect to Ed- ■ 

The next result shows that the ER reconstruction requires less sampling rate than the HR-AMP. That 
is, the ER-AMP has larger feasible operating region than the original HR-AMP under certain condition. 
Corollary 4.3: If M(ei|7?) is a concave function of ei, then we have M{dEi\r])/d < M(ei|r/). 

Proof: It is known that if a function / is concave, and /(O) > 0, then / is subadditive, i.e., 
f{x + y) < f{x) + f{y). Prom this we can get f{tx) < tf{x) for f > 1. It is clear from the definition 
that M(ei|r/) > 0. Therefore if M{Ei\r]) is concave, then by the subadditivity, we can get M(dei|??) < 
dM{Ei\r]), i.e., M{dEi\ri)/d < M{Ei\r]). ■ 

The concavity condition of M{Ei\r]) is satisfied for many families of strucfured signals. In particular, 
it is proved in 1251 for simple sparse signals in Eq. (j^ when the soft-thresholding denoiser is used. It is 
also confirmed in Q for block-sparse signals with block soft-thresholding denoiser. In the Appendix, we 
prove that it is satisfied for piecewise consfant signals. Pinally, we also show in Sec. jV] that the concavity 
condition holds for 2D images in both the simple sparse and piecewise constant families. 
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Corollary |4.3| confirms the motivation discussed in the introduction of the paper, i.e., if the CS sampling 
rate is too low, although the full-resolution reconstruction will fail, we can still reconstruct a LR version 
of the signal. Moreover, in the noiseless case, given 5i, £i, we can precisely determine the critical 
downsampling factor d hy solving the equation = M{dei\r])/d. 


C. Noise Sensitivity of MR-AMP 

The noiseless case studied above is quite restrictive. In practices, we are more interested in the 
performance of the algorithm in the presence of noise. In this part, we study the noise sensitivity of 
LR-AMP when the noises w and A(I — UdDrf)x in Eq. (j^ are not zero. As in Q, p0| , the noise 
sensitivity of HR-AMP is defined as 


Af5(c7^,(5i) = inf sup {^“(xo, (5i, r)}, 

which is the minimax MSE per coordinate of the HR-AMP output when the iteration number goes 
to cx) in Eq. ([T0|). It is shown in Q, 110| that when the undersampling ratio meets the same phase 


transition condition as in Theorem 4.1 the structured sparse signal can be recovered with a bounded 
noise sensitivity. 

When studying the noise sensitivity of the ER-AMP, we use NS{a'^^,6d) to represent the noise 
sensitivity of ER-AMP, where cr^ ^ is the variance of ER-AMP noise. The next result shows that when 
the undersampling ratio meets the same condition as in Corollary |4.2[ we can also recover the ER signal 
X(i with a bounded noise sensitivity. 

Corollary 4.4: When Cond. 2.1[ 2.2 and 2.3 are satisfied, if the undersampling ratio satisfies Eq. ([g, 
i.e., (5i > M{d£i\ri)/d in the compressed sensing of x G in Eq. ([^ with noise variance cr^, a ER 

version of the signal X(i G ^ni,ed with Ed ^ dsi can be reconstructed via ER-AMP with downsampling 
matrix and upsampling matrix U^, and the noise sensitivity is bounded by 


NS{al^,6d) 

M{d£i\r]) 


(15) 




1 - M{dEi\r])/{d6i) ^ m UrfDd)x|| 2 ). 


Proof: According to Prop. 2 in |[T0|, the noise sensitivity of AMP with various denoisers is bounded 
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by 


NSial,6i) ^ 


M{si\r]) 


1 - M{£i\r])/5i 




Replacing (5i, si and by 6d, £d and in the formula above respectively, we have 


ATQ/ 2 r N ^ M{£d\ri) ^^2 , 2 \ 

NSiad^^,6d) < ^ _ M{sd\r])/6d^ 


Since M{£d\r]) is monotonically increasing with e^, it is easy to see that i_M(£a\v)/^d 
tonically increasing. Together with < dei, we can have 


is also mono- 


NS{al^,6d) < 


M{d£i\ri) 


1 - M(dei|?7)/(d(5i) 


1^2 , 2 
[a^ + CTd,', 


By the central limit theorem, if the entries of A have i.i.d. M{0, 1/m) distribution, and and are 
deterministic, then for a given x, each entry of A(I — UrfDrf)x converges to i.i.d. Gaussian distribution 
with zero mean and variance 1/m ||(I — UdDrf)x|| 2 . Therefore the equivalent noise variance for the 
LR-AMP problem is (cj^ + 1/m ||(I — UrfDrf)x|| 2 ), which proves the result. ■ 

Different from the original AMP, the upper bound of the LR-AMP noise sensitivity NS{a‘^^,6d) is 
conditional, since it depends on the approximation error term (I — UrfDrf)x, which varies for different 
input signals. Therefore it is crucial to design good up-/down-sampling matrices to reduce the LR 
reconstruction error, which will be studied in Sec. It should be noted that the upper bound is finite 
in many applications. Moreover, sometimes we can further derive a signal-independent upper bound. For 
example, in 8-bit images, the pixel value ranges from 0 to 255. Therefore the worst value of each entry 
in (I — UrfDrf)x is 255, and the worst value of ||(I — UdDrf)x ||2 is thus 255^ni. The upper bound in 
Eq. (fTS]) can be further bounded by 


NS{al^,6d)< 

< 


M{d£i\r]) 2 255^ 

I - M{d£i\r])/{d6iy^ + 

M{d£i\r]) 2 255^(i 

1 - M{d£i\r])/{d6iy^ ^ M{d£i\r]) 


(16) 


The upper bound above is too pessimistic since the LR approximation U^Dj^x usually has much 
less approximation error than 255. The upper bound can be reduced if more accurate estimate of 
11(1 — UrfDrf)x ||2 is known. 
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Corollary |4.4| is more general than Corollary 4.2 as it allows sampling noise and LR approximation 
noise. It gives further affirmative answers to the questions raised in the introduction of the paper, i.e., if 
the CS sampling rate is too low for the full-resolution signal recovery, we can reconstruct a LR version 
of the signal with hounded noise sensitivity. The noisy case shares the same PTC with the noiseless one, 
as in the original AMP, which serves as a guideline to determine the critical resolution under which the 
noise sensitivity of the LR signal recovery is hounded. 

V. Design of Downsampling and Upsampling Matrices for MR-AMP 
In this section, we give examples on the design of the up-/down-sampling matrices that satisfy the 
three conditions in Sec. [^perfectly or approximately, so that they can he used in MR-AMP-hased image 
reconstruction. Three pairs of matrices will he designed. The first pair is in the DCT or wavelet transform 
domain and is designed for the simple sparse family. The other two pairs are in the spatial domain and 
are suitable for piecewise constant signals. 

In ITU, DCT-hased and total-variation (TV)-hased up-/down-sampling matrices are designed for videos 


such that the downsampling matrix satisfies Cond. 2.4 and the upsampling matrix satisfies Cond. 


2.1 However, fhe proof in it is mainly about TV-based up-/down-sampling matrices. Moreover, the impact 
of MR design on the quality of the measurement matrix is not considered, i.e., it is not clear whether 
Cond. 12.21 holds or not. 


A. Transform-Domain Downsampling and Upsampling 

Natural images are approximately sparse in DCT or wavelet domain. The sparse representation of a 
ni X ni image X thus belongs to the simple sparse family in Eq. and the soft-thresholding denoiser 
can be used in the transform domain. To apply CS sampling and reconstruction to images, we need to 
introduce the transform basis to Eq. Q and Eq. Q. 

Eor a ni X ni image X, a x rid ER image X^ can be obtained via transform-domain downsampling 
by first applying HR 2D transform, extracting the rid x nd low-frequency coefficients, and then applying 


the ER 2D inverse transform | |26l , 

Eet and represent the rii x ni and rid x rid DCT or orthogonal multiple-level wavelet transform 
respectively. We use the following ID transform-domain downsampling operator 
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where the fat identity matrix In^xui serves as a truncation operator, because it only keeps the first Ud 
coefficients of the input after being transformed by ■ 


Given the downsampling matrix, one way to satisfy Cond. 2.1 i.e., = I, is to use transform- 

domain zero-padding. The corresponding upsampling matrix is 


Ud = Vd^li 


m^niXUd'^nd- 


The 2D downsampling and upsampling can thus be represented as 


Xrf = DrfXDj, 
X = UrfXrfUj. 


(18) 


(19) 


It should be noted that according to the definitions in [271, for the downsampling in DCT domain, we 
can achieve non-integer downsampling ratio since we simply take the top left Ud x Ud low-frequency 
coefficients and apply the LR 2D inverse DCT transform. However, for the downsampling in wavelet 
domain, we can only get integer downsampling ratio that is power of 2, since the LR image is the 
appropriately scaled low-pass subband in the multi-level wavelet transform. 

Let X, :x.d, and x be the vectorized versions of X, X^, X, respectively, by concatenating the columns 
of each matrix together. Let (g) denote the Kronecker product, the 2D downsampling and upsampling can 
be converted to the following ID formulas. 


( 20 ) 


Xd = (Drf (g) Dd)x, 

X = (Urf (g) Urf)xrf. 

Similarly, let Si = and = IndXniSiImxnd be the 2D transform of X and its low- 

frequency part, and si and be their vectorized versions. The 2D inverse transform can be represented 
by ID transform as follows 


X = 

1 




( 21 ) 


where the two matrices are still orthogonal. Note that the corresponding ID downsampling ratio is 


n^/n^ = (P. It is easy to see that the concavity condition in Corollary 4.3 holds here, since the ID 
sparse representation of a 2D image is just the vectorized version of its 2D representation. 
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We next show that the transform-domain up-/down-sampling operators defined above satisfy Cond. 2.2 
and Cond. 12.31 


First, we assume si G ■ Since the transform-domain downsampling operator simply extracts the 

Til ,ci 

low-frequency components of si, the number of nonzero entries in is certainly no more than that in 
si; hence Sd < nfei/n^ = (fei, and E Cond. 

To check Cond. |2.2[ note that the equivalent ID measurement matrix for the HR signal is = 
(g) whereas the equivalent ID measurement matrix for the LR-CS problem in Eq. (m is 


2.3 


is thus satisfied. 




T 

nd 




( 22 ) 




xridj 


If is easy to see that ^d is the first columns of $i. Since our proposed algorithms are based on 
AMP, where each entry of the measurement matrix A follows i.i.d. AA(0, 1/m) distribution, it can be 
shown that given each entry of and ^d also has i.i.d. M{0,1/m) distribution. Therefore, with 
the proposed transform-domain up-/down-sampling method, the quality of the measurement matrix for 
the LR-AMP is the same as that of the HR-AMP. 


B. Spatial-Domain Downsampling and Upsampling 

We next develop two pairs of spatial-domain up-/down-sampling matrices for MR-AMP. In this part, 
we assume images are piecewise constant and belong to the family J^nusi Cq. which has a small 
number of change points. 

1) Solution 1: We first design the operators for ID signals and then extend them to 2D images. For 
ID piecewise constant signals, to satisfy Cond. |2.3[ the first downsampling matrix we use is the 
row-decimated identity matrix, i.e., a matrix whose {i,di)-th entries are 1 for all i, and all other entries 
are zero. The downsampled signal can be written as 


Xd = DrfX 


x[d\ x[2d] 


x[ndd] 


T 


(23) 


where x[i] represents the f-th entry of x. 


The corresponding upsampling matrix used in this part is the repetition operator which duplicates 
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each input sample by d times. 


Uh.= 


Ldxl 


Ldxl 


(24) 


where l^xi is an all-one vector. Clearly and satisfy D^Ud = I in Cond. 2.1 


Next, we show that the spatial-domain up-/down-sampling matrices also satisfy Cond. 2.3 
Lemma 5.1: If x is a piecewise constant signal generated from the family in Eq. then the 


downsampled signal in Eq. (|23| belongs to IFnFsA with Sd ^ dei. 


Ftd i^d 


Proof: 

A ni X 1 piecewise constant signal x is sparse in the differential domain. 


Si = 

-1 1 

X = ’S'„^x = 

x[2] — x[l] 

a; [3] — a; [2] 


-1 1 


- 1 

1 — 1 

1 

1 

_1 


(25) 


Similarly, the representation of the downsampling signal in the differential domain can be written as 


1 T 


(26) 


x[2d] — x[d], ..., x[ndd] — x[{nd — l)d] 

Therefore, calculating the number of change points in x^ is equivalent to counting the number of nonzero 
entries in s^. 


x[l] sf 


and s* = 


1 T 


x[d] s 


, i.e., 


To facilitate the proof, we construct two new vectors = 
adding the first entry of x and x^ to si and respectively. If we add d consecutive entries of s^, we can get 
one entry of s^. Eor example, {x[d+l]—x[d]) + {x[d+2]—x[d+l]) + ... + {x[2d]—x[2d—l]) = x[2d]—x[d]. 
In matrix form, this means = Ujs^. If x is generated from the maximum expected number 

of nonzero entries in will be ni£i -|- 1, due to the extra x[l] in it. According to = Ujs^, the 
maximum expected number of nonzero entries in is still niei -|- 1. This happens when there is at 
most one nonzero entry in every d entries in s^; hence £d < (niei -\-l)/nd = dei + l/nd —)• dsi when 


ni —>■ oo. 
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We next extend the results above to 2D images. The 2D Ud x Ud LR image can be written as Eq. 
( [T^ with Dd in Eq. ( |2^ . If an image is piecewise constant, its 2D gradient is sparse, where the 2D 
gradient at each pixel is given by 


(VX)jj — Xjj]. 


(27) 


The number of change points in a 2D piecewise constant signal X equals to the number of nonzero 
entries in VX, where (VX)jj is counted as one nonzero entry if one or two of its components are 
nonzero. Therefore, we can also vectorize the 2D VX into a ID vector, and apply the method in the 


Appendix to prove the concavity in Corollary 4.3 for 2D piecewise constant signals. Additionally, the 


vertical differences and the horizontal differences are disjoint. By Eemma [53 the number of horizontal 
or vertical change points of X^ is no larger than that of X, thus Cond. |2.3| is true for 2D images. 

The remaining problem is to choose the appropriate denoiser for 2D piecewise constant signals. In this 


paper, instead of using the denoisers discussed in 1101, |211, such as NEM (non-local means) and BM3D 
(3D block matching), we use a 2D-TV-based denoiser in r/o-*(z^) of Eq. (j^. Our method is denoted as 
AMP-TV-2D. 

The TV norm of 2D piecewise constant signals is defined as 


^IItV - ^ y v + ’ 


(28) 




which is isotropic and un-differentiable. This norm will be used by the 2D-TV-based denoiser. Eurther 


details are given in Sec. VI-A It is different from the ID TV denoiser in Q where TV norm for ID 

ni —1 

piecewise constant signal is written as ||x||.py = ^ — Xi\. 


i=l 


In Sec. VI-D4[ we compare the performance of our AMP-TV-2D with the state-of-the-art algorithm 
TVAE3 (TV minimization by Augmented Eagrangian and AEternating direction AEgorithms) in | |28| . 
Note that TVAE3 depends on two slack parameters, which have to be manually tuned for each image and 
each measurement rate. In contrast, the thresholding parameters in our AMP-TV-2D are automatically 
tuned in each iteration, which will be discussed in Sec. VI-A| Recently, a similar algorithm to our AMP- 


TV-2D, the dual-constraints AMP (DC-AMP) (Sec. 8.1 of |24|), is developed for 2D piecewise smooth 
signals, which can achieve similar performance to TVAE3. However, it also has a smoothness parameter 
that needs to be manually tuned. Moreover, there is no theoretical analysis for DC-AMP. 
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Given the spatial-domain up-/down-sampling matrices, to satisfy Cond. 2.2 i.e., the quality of the mea¬ 
surement matrix for LR-AMP is no worse than that of HR-AMP, we need to normalize the measurement 
matrix for LR-AMP, i.e., 

^ (29) 


= ^A(Ud <8) Ud), 

such that each entry of has i.i.d. Af{0,1/m) distrihution. 

2) Solution 2: In addition to the simple up-/down-sampling matrices in Eq. 


and Eq. (24i, we 


also develop a pair of hicuhic upAdownsampling matrices and evaluate them in Sec. VI-D In hicuhic 
downsampling, each pixel in the ER image is the weighted average of sixteen pixels in the HR image. 


which has been known to produce smoother ER image than Eq. (231, i.e., with less number of change 


points in X^. Therefore Cond. |23] holds for bicubic downsampling. On the other hand, the upsampling first 
inserts d—1 zeros between neighboring samples of the ER image and then performs bicubic interpolation. 
However, it can be verified that the corresponding product D^Ud is not an identity matrix, although very 


close. Therefore, strictly speaking, Cond. 2.2 does not hold for bicubic matrices, and the simple scaling 
matrix A cannot make each entry of exactly having i.i.d. A((0,1/m) distribution. Nevertheless, this 


is still approximately true, and the efficiency of this scheme will be verified empirically in Sec. VI-D 


Moreover, according to Corollary |4.4[ the conditional upper bound of the noise sensitivity is proportional 
to the ER approximation error ||(I — UrfDrf)x|| 2 . Therefore for images, in terms of ER approximation 


error, the bicubic up-/down-sampling matrices are still better than the simple matrices in Eq. (231 and 
Eq. 

Einally, we point out the differences of our methods with those in 1141, |16|. In | [T4] |, a similar spatial- 
domain up-/down-sampling framework was proposed, but the proof in it was implicit. Also, TVAE3 
was chosen as the reconstruction algorithm, which requires manual tuning of two parameters. Moreover, 
the reconstruction performance cannot be predicted. Our AMP-TV-2D does not have manually tuned 


parameter, and its performance can be accurately predicted via state evolution. In 116|, the same piecewise 
constancy model and the up-/down-sampling matrices as in Eq. ( [23] ) and (24i are used. It first reconstructs 
the original HR image and uses this estimated HR image to reduce the approximation error A(I—UrfDrf)x. 
However, there is no theoretical guarantee that such operation can reduce the approximation error, and 
the algorithm only works when the undersampling rate is sufficiently large, at least 20%. Moreover, 
the complexity of this approach is higher than reconstructing the ER image directly. 
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Fig. 1. Empirical intermediate MSE and predicted state evolution of HR-AMP-ST and LR-AMP-ST for image Barbara with 

d = 4. 


VI. Experimental Results 


In this section, we demonstrate the performance of the proposed MR-AMP with both transform-Zspatial- 
domain up-/down-sampIing, denoted hy AMP-ST (soft thresholding) and AMP-TV (total variation), 
respectively. Empirical results will also he shown to verify some theoretical results. In each method, to 
facilitate comparison with the conventional approach, we use ER-AMP-ST and ER-AMP-TV to denote 
the proposed ER reconstruction schemes, and HR-AMP-ST and HR-AMP-TV to denote the original AMP 
with HR reconstruction. In addition, H2E-AMP-ST and H2E-AMP-TV represent the naive solutions that 
first reconstruct the HR signal and then downsample to the ER. 

All tests in this paper use column-normalized i.i.d. Gaussian measurement matrix A. All simulations 
are conducted on a PC with 3.4GHz Intel Core i7 quad-core processor and 64GB of memory. The testing 
images used include popular images Eena, Barbara, Boat, House, and Peppers, as well as some land 
remote sensing images, including the Memorial Stadium at the University of Nebraska Cornhuskers, and 


Sea World in San Diego. We follow the setup in 110| to rescale all images to 128 x 128. This enables the 
entire measurement matrix A to be stored in the memory. We also include some experiments of larger 


256 X 256 images to demonstrate the visual comparison, following the same setup in |10|. 


A. Parameter Tuning 

One of the main challenges in implementing different MR-AMP algorithms is the tuning of each 
algorithm’s free parameters. Many techniques exist to estimate the noise variance in an image. In this 
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• Predicted HR-AMP-TV-2D, d=1 

■ Observed HR-AMP-TV-2D, d=1 

• Predicted LR-AMP-TV-2D, d=2 

■ Observed LR-AMP-TV-2D, d=2 

• Predicted LR-AMP-TV-2D, d=4 

■ Observed LR-AMP-TV-2D, d=4 




••r‘rr**rr a 


15 

Iteration 



-Predicted HR-AMP-TV-2D, d=1 

-Observed HR-AMP-TV-2D, d=1 

- Predicted LR-AMP-TV-2D, d=2 

- Observed LR-AMP-TV-2D, d=2 

- Predicted LR-AMP-TV-2D, d=4 

- Observed LR-AMP-TV-2D, d=4 




15 

iteration 


(a) (b) 

Fig. 2. State evolutions of MR-AMP-TV with a CS sampling rate of 5% and no measurement noise for the 128 x 128 Barbara 
image, (a) Repetition interpolator, (b) Bicubic interpolator. 


paper, we use the following convenient feature of AMP algorithms: llr^ll^/mRi (cr^)^ i29|. 

For MR-AMP-ST, we set its threshold using three methods. For the ID synthetic examples in Sec. 
VI-C[ we assume the sparsity rate is known and set the thresholding parameter according to the minimax 


rule in | [TT[ . For the 2D imaging examples in Sec. VI-D[ since images are not exactly sparse in transform 
domain, we have to estimate the sparsity rate. For sufficient large CS undersampling rate such as 
10% and 20%, we use the SURE (Stein’s unbiased risk estimate)-based method in | [^ to decide the 
thresholding parameter in each iteration. For very small (5i such as 3% and 4%, SURE does not work 
well since it is based on large system limit, we choose the max-min optimal threshold as determined by 

ED- 

For AMP-TV, we use different tuning methods for ID and 2D signals. For ID signals, we use the 
source code from |[3^ directly. For 2D images, there are many methods on how to adaptively choose 


the regularization parameter in TV-based image denoising, e.g., |33| and |34|. In this paper, we use 


Algorithm 6 in |34|, due to its simplicity and efficiency. In each iteration of AMP-TV-2D, a Lagrangian 
optimization problem is solved, whose constraint is the TV of the solution, and the Lagrangian parameter 
can be adaptively determined by a formula. 

In AMP-ST, the Onsager term is obtained by Eq. (4.1) in Q. For AMP-TV-ID, the Onsager term 
is calculated by Eq. (5.11) in Q. Eor AMP-TV-2D, it is difficult to obtain an exact expression of the 


divergence. We thus apply the Monte Carlo method in 1101 to hnd a good approximation of the divergence. 
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B. State Evolution in MR-AMP 

In this part, we compare the predicted and observed performances of MR-AMP with different denoisers. 
Recall that the state evolution of MR-AMP is given in Eq. To compute this value, at every iteration 
we add white Gaussian noise with standard deviation cr^ to o> denoise the signal with denoiser (:, r), 
and then compute the MSE. 

Eig. [T] compares the empirical MSE and predicted state evolution of MR-AMP-ST for the test image 
Barbara of size 128 x 128, with DCT being the sparsifying basis. It can be seen that the state evolution 
is quite accurate. Moreover, the converged MSE per entry of the ER image is about 50% smaller than 
that of the HR image, which verifies the motivation of this paper, i.e., we can recover a ER signal with 
smaller MSE when the MSE of the HR signal is too large. Note that the ER reference image is obtained 
via the DCT-domain downsampling in Sec. V-A[ and the corresponding MSE is the MSE between the 
reconstructed ER image by ER-AMP-ST and the ER reference image. 

Eig. 1^ shows the state evolution performance of MR-AMP-TV. Two different upsampling matrices 


are compared: the repetition interpolator in Eq. (24i (MR-AMP-TV-2D-R) and bicubic interpolator 
(MR-AMP-TV-2D-B). The reference ER image is obtained by Matlab’s imresize(x,l/d) command with 
bicubic interpolator. There is near perfect correspondence between the predicted and true MSEs for the 
repetition interpolation. Eor bicubic interpolator, a slight mismatch exists, because the entries of the new 
measurement matrix are not exactly independent. The figures also show fhaf lower resolution provides 
smaller MSE, and bicubic interpolator outperforms the repetition operator. 

Note that the denoiser in the AMP-TV-2D is essentially a non-scalar denoiser, similar to p^ , |j^, 
1|24). Although the state evolution for AMP with non-scalar denoisers has not been proved rigorously, 
the results in Pig. suggest that the state evolution derived in our paper is quite accurate. 


C. Performance with Synthetic ID Signals 

In this part, we demonstrate the performance of the proposed scheme for synthetic ID signals, which 
can verify the theoretical noiseless phase transition curve (PTC) and noise sensitivity. 

1) Transform Domain Approach: To get the empirical noiseless PTC of HR-AMP-ST, we fix rei = 
2000, and take 30 equally distanced values of (5i = m/m in the range of [0.05,0.95], and 30 equally 
distanced values of pi = ki/m in [0.05,0.95]. Eor each combination of (5i, pi), a ID Bernoulli-Gaussian 
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7 

HR-AMP-ST 

Bound 

HR-AMP-ST 

Empirical 

ER-AMP-ST 

Bound 

ER-AMP-ST 

Empirical 

0.95 

3.80 

3.06 

5.23 

2.78 

0.98 

9.80 

7.47 

5.23 

4.12 

0.99 

19.80 

14.62 

5.23 

4.15 

0.998 

39.80 

28.89 

5.23 

4.79 


TABLE I 

Noise sensitivity of MR-AMP-ST with 5i = 0.2 and pi = 0.3. 


signal and its CS samples are generated before applying the HR-AMP-ST. The empirical PTCs are 
obtained by connecting operating points with 50% success rate of the signal recovery, where the recovery 

Il2/ll Il2 

is considered successful when the normalized MSE (NMSE) satisfies ||xo — x||2 / ||xo II2 ^ 10" ■ 

To study the empirical noiseless PTC of ER-AMP-ST, we generate a special ni x 1 sparse signal, 
whose first Ud = ni/d entries are Bernoulli-Gaussian distributed, and all other entries are 0. According 
to Eq. ([^, the truncation operator does not introduce any approximation error A(I — UrfDrf)x. We then 
run the HR-AMP-ST and ER-AMP-ST algorithms to recover the target HR and ER signals respectively. 
Note that we are interested in the case m/ni < 1/d, otherwise the setup is no longer a CS problem. 
Although the procedure of generating the HR signal here is different from that in the simulation of 
empirical HR-AMP-ST above, both signals belong to the same class of probability distribution if the 
numbers of nonzero coefficients are the same, and experimental results show that these two empirical 
PTCs for HR-AMP-ST coincide with each other. 

The theoretical noiseless PTC in Eq. ( [14] ) and the empirical noiseless PTC of ER-AMP-ST are shown 
in Pig. for simple sparse signals with different d. The two sets of curves agree perfectly. It can be 
shown that as d increases, the PTC curve shifts to the left, which means that the ER-AMP can recover 
the signal even when the HR-AMP fails. 

The example above does not have approximation error. Next, we construct a special case to show that 
the noise sensitivity of HR-AMP-ST is unbounded above the PTC, while the noise sensitivity of the 
ER-AMP-ST is still bounded. The setup is similar to that in Q, where a special 3-point distribution of x 
is constructed in Eemma 4.4, whose MSE above the phase transition boundary is given by < 5 i 7 /(l — 7 ). 
Therefore the MSE can go to infinity when 7 is close to 1. We present in Table [^ the noisy sensitivity of 
MR-AMP-ST with ni = 2000, (5i = 0.2, pi = 0.3 and = 1. As shown in Pig. this setup is above 
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Fig. 3. The theoretical and empirical PTCs of MR-AMP-ST. 


the PTC of d = 1, but below the PTC of d = 2. The non-zero locations of x are chosen with probability 
1.8ei from the first n 2 entries to generate the 3-point distribution, and with probability 0.2ei to generate 
Bemoulli-Gaussian signals for the second n 2 entries, in order to fix the approximation error in Eq. Q 
for different 7 ’s. We then apply HR-AMP-ST and LR-AMP-ST to reconstruct x and x^. 

It can be seen from Table that as 7 approaches to 1, the noise sensitivity bound of HR-AMP-ST 


keeps increasing, but the noise sensitivity bound of LR-AMP-ST is stable because all parts in Eq. (151 
are fixed. This verifies the advantage of our LR-AMP. The empirical results of both methods are also 
below their noise sensitivity bounds. 

2) Spatial Domain Approach: It is difficult to reproduce the theoretical noiseless PTC of HR-AMP- 
TV-ID in Q since it relies on complicated numerical optimization and no open source code is available. 
Instead, we study the empirical noiseless PTC of HR-AMP-TV-ID by replicating an experiment from 


1321 using its source code. We fix ni = 628, and consider a 30 x 30 uniform grid in the range of 
= m/ni G [0.05,0.95] and pi = ki/m G [0.05,0.95]. The corresponding HR Bernoulli-Gaussian ID 
finite-difference signal is then generated. The empirical noiseless PTC of HR-AMP-TV-ID is shown in 
Pig. (a) (with d = l). 

To get the empirical noiseless PTCs of LR-AMP-TV-ID, we first generate LR signal x^ that yields 
ID Bernoulli-Gaussian finite-difference sequence with sparsity rate dei. We then duplicate each entry d 


times to get the HR piecewise constant signal with sparsity rate ei, according to and in Eq. (231 


and (24i. Prom the analysis in Sec. V-B[ the approximation error is zero. Successful recovery is declared 
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(a) (b) 


Fig. 4. (a) The empirical PTCs of MR-AMP-TV-ID for Bernoulli-Gaussian finite-difference signals, (b) MR recovery of 

Bemoulli-Gaussian finite-difference signals with sparsity rate e\ = 0.05 and SNR of 60dB in the measurement. 



(a) (b) (c) (d) (e) (f) 

Fig. 5. Reconstruction of 10% sampled 256 x 256 Barbara image with dowansampling factor d = 2 and DCT as the 
sparsifying basis for MR-AMP-ST. (a) HR-AMP-ST (20.32dB). (b) H2L-AMP-ST (21.31dB). (c) LR-AMP-ST (22.72dB). (d) 
HR-AMP-TV-2D (25.06dB). (e) H2L-AMP-TV-2D (27.75dB). (f) LR-AMP-TV-2D-B (27.54dB). 


when NMSE is below 10“^. The results with d = 2 and d = 4 are also shown in Fig. |^(a). 

To study the noise sensitivity of MR-AMP-TV-ID, we recover the target HR and LR piecewise constant 
signals after introducing additional white Gaussian noise (AWGN) with SNR = || Ax ||2 /||w ||2 = 60dB 
in the measurement. Fig. [^(b) shows the median NSNR defined as NSNR = ||xo ||2 / ||xo — x ||2 versus 
sampling ratio <5i = m/ni at the fixed sparsity rate ei = 0.05, as in 1^. It shows that the FR-AMP-TV- 


ID has lower NMSE than the HR-AMP-TV-ID. This verifies Corollary |4.4[ i.e., the FR reconstruction 
has better performance than the HR one. 


D. Performance with 2D Images 

In this part, we apply the MR-AMP theory to MR 2D image reconstruction. All reported experimental 
results are the averages of 20 Monte Carlo simulations. 
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d 


Algorithm 

Eena 

Barbara 

Boat 

House 

Peppers 

HuskerStadium 

SeaWorld 



HR-AMP-ST 

16.75 

15.96 

17.60 

18.21 

15.53 

15.86 

14.39 


5% 

H2E-AMP-ST 

17.40 

16.48 

18.30 

18.58 

15.93 

16.49 

15.10 



ER-AMP-ST 

18.02 

17.11 

18.77 

19.13 

16.68 

16.89 

15.33 



HR-AMP-ST 

18.50 

17.79 

18.94 

19.71 

17.35 

16.95 

15.17 


10% 

H2E-AMP-ST 

19.43 

18.56 

19.97 

20.34 

18.19 

18.05 

16.10 

9 


ER-AMP-ST 

20.82 

19.94 

21.07 

21.72 

19.43 

18.71 

16.79 



HR-AMP-ST 

21.28 

20.36 

21.08 

22.31 

19.93 

18.69 

16.60 


20% 

H2E-AMP-ST 

22.58 

21.69 

22.61 

23.46 

21.27 

20.13 

18.06 



ER-AMP-ST 

24.90 

24.25 

24.46 

26.34 

23.76 

21.89 

19.72 



HR-AMP-ST 

15.37 

14.76 

16.54 

17.09 

14.33 

14.96 

13.51 


3% 

H2E-AMP-ST 

16.98 

16.33 

18.55 

18.86 

15.91 

17.24 

15.97 



ER-AMP-ST 

18.22 

17.66 

18.92 

19.58 

17.02 

17.13 

15.59 



HR-AMP-ST 

16.03 

15.24 

16.95 

17.56 

14.87 

15.27 

13.75 


4% 

H2E-AMP-ST 

17.81 

16.91 

19.09 

19.46 

16.65 

17.71 

16.39 

4 


ER-AMP-ST 

19.21 

18.42 

19.63 

20.31 

17.78 

17.66 

15.76 


HR-AMP-ST 

16.52 

15.74 

17.24 

17.97 

15.91 

15.52 

13.95 


5% 

H2E-AMP-ST 

18.44 

17.54 

19.57 

20.04 

17.29 

18.13 

16.72 



ER-AMP-ST 

19.66 

18.90 

19.67 

20.60 

18.45 

17.48 

15.72 


TABLE II 

PSNRS (DB) of 128 X 128 IMAGE RECONSTRUCTIONS WITH DCT-DOMAIN MR-AMP-ST. 






(a) (b) (c) (d) (e) 

Fig. 6. Reconstructed of 20% sampled 256 x 256 HuskerStadium image with downsampling factor d = 2 and D8 wavelet as 
the sparsifying basis for MR-AMP-ST. (a) HR-AMP-ST (18.65dB). (b) H2L-AMP-ST (20.22dB). (c) LR-AMP-ST (21.05dB). 
(d) HR-AMP-TV-2D (21.66dB). (e) H2L-AMP-TV-2D (24.52dB). (f) LR-AMP-TV-2D-B (24.38dB) . 


1) Target LR image: The target LR images are different when different downsampling matrices are 
used. For the transform-domain approach, the target LR image is represented by Eq. ( [T^ . Both DCT 
and the Daubechies-8 (D8) wavelet are tested. For the spatial-domain approach, although the simple 


matrix in Eq. (231 can be applied, we choose to use the bicubic downsampling matrix, as it leads to 


better ER image. As discussed before, Cond. 2.3 still holds in this case. Given the bicubic downsampling 
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d 


Algorithm 

Eena 

Barbara 

Boat 

House 

Peppers 

HuskerStadium 

SeaWorld 



HR-AMP-ST 

16.58 

15.84 

17.81 

17.66 

15.31 

15.94 

14.33 


5% 

H2E-AMP-ST 

16.85 

16.46 

18.55 

18.30 

15.83 

16.86 

15.04 



ER-AMP-ST 

17.35 

17.01 

19.13 

18.95 

16.67 

17.28 

15.38 



HR-AMP-ST 

18.20 

17.47 

19.29 

19.56 

17.14 

16.99 

15.05 


10% 

H2E-AMP-ST 

19.13 

18.25 

20.53 

20.43 

18.05 

18.18 

15.98 

9 


ER-AMP-ST 

20.62 

19.72 

21.46 

21.97 

19.46 

19.02 

16.79 



HR-AMP-ST 

21.27 

20.15 

21.62 

22.68 

20.04 

18.86 

16.60 


20% 

H2E-AMP-ST 

22.98 

21.59 

23.53 

24.47 

21.61 

20.64 

18.11 



ER-AMP-ST 

24.98 

23.89 

25.02 

26.44 

23.51 

22.05 

19.80 



HR-AMP-ST 

15.14 

14.67 

16.66 

16.41 

14.26 

15.01 

13.54 


3% 

H2E-AMP-ST 

16.83 

16.40 

18.90 

18.31 

15.95 

17.45 

16.19 



ER-AMP-ST 

17.83 

17.24 

19.31 

19.28 

16.93 

17.33 

15.52 



HR-AMP-ST 

15.62 

15.16 

17.09 

17.08 

14.70 

15.40 

13.69 


4% 

H2E-AMP-ST 

17.47 

17.04 

19.53 

19.21 

16.63 

18.07 

16.50 

4 


ER-AMP-ST 

18.73 

18.09 

19.70 

19.78 

17.60 

17.74 

15.78 


HR-AMP-ST 

15.99 

15.58 

17.50 

17.52 

15.19 

15.66 

13.90 


5% 

H2E-AMP-ST 

17.99 

17.63 

20.18 

19.81 

17.25 

18.51 

16.87 



ER-AMP-ST 

19.00 

18.47 

19.65 

20.14 

17.84 

17.55 

15.60 


TABLE III 

PSNRS (DB) of 128 X 128 IMAGE RECONSTRUCTIONS WITH WAVELET-DOMAIN MR-AMP-ST. 


matrix, we test the repetition upsampling matrix in Eq. ( [24| ) as well as the hicuhic upsampling matrix. It 
can he verified that Cond. |2.1| = I holds approximately between these two upsampling matrices 

and the hicuhic downsampling matrix. 

In this paper, we use the Peak SNR (PSNR) to measure the objective quality of a reconstructed image, 
which is defined as 101og2o(255^/MSE(X — X)), where X is the reference image, and X is the test 
image. 


2) Scaling Matrix A: During the reconstruction of LR image, in order to ensure that Cond. 2.2 


m 


Sec. 1^ is satisfied, we need to scale its corresponding measurement matrix AU^ into = AU^A 
to get normalized columns, as shown in Eq. ( [22| ) and Eq. pQ] ). Since no specific entries in the target 
ER image are preferred, the scaling matrix A should be diagonal matrix with equal diagonal entries. 
Eor ER-AMP-ST in DCT and wavelet domain, the diagonal entry is the inverse of the downsampling 
factor d, according to Eq. ( [22] ). Eor ER-AMP-TV-2D in TV domain, things are slightly different. Eor the 
repetition operator that replaces each pixel in ER image with a. d x d block of pixels in the HR image. 
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d 

Si 

Algorithm 

Lena 

Barbara 

Boat 

House 

Peppers 

HuskerStadium 

SeaWorld 



HR-AMP-TV-2D 

20.88 

19.66 

20.67 

22.85 

19.60 

18.36 

16.13 


5% 

H2L-AMP-TV-2D 

22.55 

21.09 

22.51 

24.53 

21.07 

20.27 

17.93 


LR-AMP-TV-2D-R 

21.79 

20.25 

22.17 

23.80 

20.13 

19.95 

17.75 



LR-AMP-TV-2D-B 

22.68 

21.17 

22.67 

25.07 

21.13 

20.33 

17.87 



HR-AMP-TV-2D 

23.62 

22.18 

22.80 

26.55 

22.51 

20.20 

17.64 


10% 

H2L-AMP-TV-2D 

25.83 

24.08 

25.35 

28.91 

24.63 

22.70 

19.97 

z 

LR-AMP-TV-2D-R 

23.83 

22.32 

24.07 

26.42 

22.31 

21.64 

19.16 



LR-AMP-TV-2D-B 

25.66 

24.18 

25.24 

28.82 

24.33 

22.63 

19.95 



HR-AMP-TV-2D 

26.51 

25.05 

25.02 

30.91 

25.70 

22.15 

19.31 


20% 

H2L-AMP-TV-2D 

29.49 

27.65 

28.39 

34.11 

28.49 

25.38 

22.26 


LR-AMP-TV-2D-R 

26.24 

24.83 

26.19 

29.23 

24.60 

23.60 

20.85 



LR-AMP-TV-2D-B 

28.92 

27.63 

27.99 

32.44 

27.51 

25.28 

22.34 



HR-AMP-TV-2D 

18.69 

17.90 

18.89 

20.32 

17.71 

16.77 

15.08 


3% 

H2L-AMP-TV-2D 

21.75 

20.80 

22.29 

23.55 

20.70 

20.26 

18.71 


LR-AMP-TV-2D-R 

20.73 

19.63 

22.10 

23.59 

19.33 

20.22 

18.80 



LR-AMP-TV-2D-B 

21.95 

20.49 

23.02 

24.47 

20.46 

20.92 

19.22 



HR-AMP-TV-2D 

19.89 

18.89 

19.84 

21.64 

18.83 

17.66 

15.65 

A 

4% 

H2L-AMP-TV-2D 

23.43 

22.21 

23.77 

25.39 

22.26 

21.62 

19.75 


LR-AMP-TV-2D-R 

21.53 

20.28 

22.75 

24.19 

20.31 

20.81 

19.35 



LR-AMP-TV-2D-B 

22.99 

21.54 

23.93 

25.55 

21.66 

21.64 

19.93 



HR-AMP-TV-2D 

20.88 

19.66 

20.67 

22.85 

19.60 

18.36 

16.13 


5% 

H2L-AMP-TV-2D 

24.86 

23.24 

25.07 

27.05 

23.40 

22.74 

20.49 


LR-AMP-TV-2D-R 

22.24 

20.86 

23.41 

25.09 

20.75 

21.51 

19.66 



LR-AMP-TV-2D-B 

23.97 

22.32 

24.80 

26.58 

22.52 

22.54 

20.30 


TABLE IV 

PSNRS (DB) of 128 X 128 IMAGE RECONSTRUCTIONS WITH SPATIAL-DOMAIN MR-AMP-TV-2D. 


the diagonal entry in the scaling matrix is still 1/d. For hicuhic interpolation, we empirically set the 
diagonal entry in the scaling matrix to he 1/2.68 for d = 2 and 1/5 for d = 4. Although this approach 
cannot exactly normalize the columns and there are still some correlations between entries in the new 
measurement matrix, it works quite well in practices. 

3) Noiseless image recovery: Tables |T^ |T^ and IV compare the performances of DCT-domain MR- 
AMP-ST, wavelet-domain MR-AMP-ST, and spatial-domain MR-AMP-TV when there is no measurement 
noise. In each case, we compare our proposed LR-AMP that recovers the LR image directly, the 
conventional HR-AMP that reconstructs the HR image, and the naive H2L-AMP that recovers the HR 
image first and then downsamples it to obtain the LR image with the corresponding downsampling matrix. 
The highest PSNR in each case is highlighted. 
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d 


Algorithm 

Lena 

Barbara 

Boat 

House 

Peppers 

HuskerStadium 

SeaWorld 

2 

10% 

HR-AMP-TV-2D 

HR-TVAL3 

23.62 

23.36 

22.18 

21.80 

22.80 

22.94 

26.55 

26.21 

22.51 

21.89 

20.20 

20.32 

17.64 

17.71 

LR-AMP-TV-2D-B 

LR-TVAL3 

25.66 

25.44 

24.18 

24.05 

25.24 

25.09 

28.82 

28.60 

24.33 

23.98 

22.63 

21.94 

19.95 

19.24 

20% 

HR-AMP-TV-2D 

HR-TVAL3 

26.51 

26.80 

25.05 

25.22 

25.02 

25.49 

30.91 

31.79 

25.70 

25.65 

22.15 

22.42 

19.31 

19.62 

LR-AMP-TV-2D-B 

LR-TVAL3 

28.92 

28.58 

27.63 

27.45 

27.99 

27.59 

32.44 

31.67 

27.51 

27.02 

25.28 

24.19 

22.34 

21.38 

4 

4% 

HR-AMP-TV-2D 

HR-TVAL3 

19.89 

19.69 

18.89 

18.77 

19.84 

20.21 

21.64 

21.43 

18.83 

18.48 

17.66 

18.21 

15.65 

16.02 

LR-AMP-TV-2D-B 

LR-TVAL3 

22.99 

22.77 

21.54 

21.31 

23.93 

23.85 

25.55 

25.19 

21.66 

21.56 

21.64 

20.95 

19.93 

19.80 

5% 

HR-AMP-TV-2D 

HR-TVAL3 

20.88 

20.58 

19.66 

19.29 

20.67 

20.94 

22.85 

22.57 

19.60 

19.15 

18.36 

18.70 

16.13 

16.38 

LR-AMP-TV-2D-B 

LR-TVAL3 

23.97 

23.80 

22.32 

21.85 

24.80 

24.40 

26.58 

26.15 

22.52 

22.46 

22.54 

21.45 

20.30 

20.45 


TABLE V 

Comparison of the final reconstruction results in PSNR between TVAL3 and AMP-TV-2D. 


From Tables and III we can see that LR-AMP-ST almost always outperforms the other two 
algorithms, except when d = 4 for HuskerStadium and SeaWorld. This is partially due to two reasons. 
First, land remote sensing images contain more details compared to natural images. Second, the suboptimal 


thresholding rule in [311 is used for d = A, whereas the optimal SURE-based thresholding method in 
| [^ is used for d = 2. 

In the spatial-domain approach, LR-AMP-TV-2D-B and H2L-AMP-TV-2D are the top two algorithms. 
Their reconstruction performances are comparable and the PSNR difference between them is within 1 
dB. However, H2L-AMP-TV-2D is much slower than the proposed LR-AMP-TV-2D, as detailed in the 
computational complexity part later. Since the reference HR image is the same for the three approaches 
listed in Tables and IV it can be seen that the TV-based approach yields higher PSNR than the 

transform-domain ones. 

Fig. and Fig. illustrate the visual quality of the recovered 256 x 256 Barbara and Stadium by 


different methods. It can be seen that transform-domain and spatial-domain approaches have different 
types of reconstruction artifacts. The former preserves more details but also contains more high frequency 
noises, whereas the latter is blockier, despite higher PSNRs. 
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AWGN with standard deviation 20 





5% 

10% 

15% 


<^1 

3% 

4% 

5% 

DCT 

r\ 1 

HR-AMP-ST 

16.00 

17.70 

19.92 

^ A 

HR-AMP-ST 

14.74 

15.22 

15.73 

Q — Z 

H2L-AMP-ST 

16.45 

18.43 

21.08 

Q—4 

H2L-AMP-ST 

16.32 

16.90 

17.53 



LR-AMP-ST 

17.12 

19.65 

22.88 


LR-AMP-ST 

17.56 

18.37 

18.71 




5% 

10% 

15% 


<5i 

3% 

4% 

5% 

Wavelet 


HR-AMP-ST 

15.80 

17.47 

19.64 


HR-AMP-ST 

14.65 

15.14 

15.56 

d=2 

H2L-AMP-ST 

16.44 

18.33 

21.02 

d=4 

H2L-AMP-ST 

16.38 

17.04 

17.61 



LR-AMP-ST 

16.85 

19.56 

22.79 


LR-AMP-ST 

17.20 

17.98 

18.29 




5% 

10% 

15% 


<^1 

3% 

4% 

5% 



HR-AMP-TV-2D 

19.59 

21.84 

23.93 


HR-AMP-TV-2D 

17.82 

18.85 

19.59 

TV 

d=2 

H2L-AMP-TV-2D 

21.00 

23.68 

26.24 

d=4 

H2L-AMP-TV-2D 

20.72 

22.13 

23.09 



LR-AMP-TV-2D-R 

20.20 

22.08 

24.17 


LR-AMP-TV-2D-R 

19.60 

20.25 

20.80 



LR-AMP-TV-2D-B 

21.05 

23.70 

26.31 


LR-AMP-TV-2D-B 

20.45 

21.48 

22.22 


AWGN with standard deviation 40 





5% 

10% 

15% 


<^1 

3% 

4% 

5% 

DCT 

r\ 1 

HR-AMP-ST 

15.89 

17.37 

19.06 

A A 

HR-AMP-ST 

14.67 

15.16 

15.59 

Q — Z 

H2L-AMP-ST 

16.34 

18.03 

19.94 


H2L-AMP-ST 

16.26 

16.96 

17.43 



LR-AMP-ST 

17.00 

19.11 

21.21 


LR-AMP-ST 

17.41 

18.09 

18.24 



<5i 

5% 

10% 

15% 


<^1 

3% 

4% 

5% 

Wavelet 

A 1 

HR-AMP-ST 

15.68 

17.19 

18.72 

A A 

HR-AMP-ST 

14.61 

15.08 

15.42 

U—Z 

H2L-AMP-ST 

16.21 

17.87 

19.75 

Q—4 

H2L-AMP-ST 

16.37 

17.00 

17.49 



LR-AMP-ST 

16.81 

18.95 

21.10 


LR-AMP-ST 

17.00 

17.76 

17.93 




5% 

10% 

15% 


<5i 

3% 

4% 

5% 



HR-AMP-TV-2D 

19.36 

21.13 

22.52 


HR-AMP-TV-2D 

17.75 

18.70 

19.36 

TV 

d=2 

H2L-AMP-TV-2D 

20.70 

22.80 

24.51 

d=4 

H2L-AMP-TV-2D 

20.62 

21.94 

22.75 



LR-AMP-TV-2D-R 

20.02 

21.58 

23.05 


LR-AMP-TV-2D-R 

19.55 

20.16 

20.61 



LR-AMP-TV-2D-B 

20.73 

22.81 

24.44 


LR-AMP-TV-2D-B 

20.35 

21.31 

21.91 


TABLE VI 

PSNRs (dB) of reconstruction of 128 X 128 Barbara image with varying amounts of additive Gaussian 

MEASUREMENT NOISE. 


4) Comparison between AMP-TV-2D-B and optimal TVAL3: In Table |Vj we compare the results of 
TVAL3 with optimized slack parameters |[T4|, |[2^ and our parameter-free AMP-TV-2D-B for the MR-CS 


problem in Eq. Q- For the original HR image reconstruction, the performance of HR-AMP-TV-2D-B 
is comparable to the optimized HR-TVAL3. However, for the LR image reconstruction, our LR-AMP- 
TV-2D outperforms the optimized LR-TVAL3 in almost all cases by up to IdB. More importantly, the 


theoretical analyses developed in Sec. IV and |V] are applicable for MR-AMP-TV, whereas there are only 
some qualitative analyses in |[T4t. 
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d 


Algorithm 

Lena 

Barbara 

Boat 

House 

Peppers 

HuskerStadium 

SeaWorld 



HR-AMP-ST 

18.50 

17.79 

18.94 

19.71 

17.35 

16.95 

15.17 

9 

10% 

L2H-AMP-ST 

20.14 

19.46 

20.18 

21.25 

19.02 

18.03 

16.07 


HR-AMP-TV-2D-B 

23.62 

22.18 

22.80 

26.55 

22.51 

20.20 

17.64 



L2H-AMP-TV-2D-B 

23.65 

22.52 

22.86 

26.21 

22.30 

20.25 

17.79 



HR-AMP-ST 

16.52 

15.74 

17.24 

17.97 

15.91 

15.52 

13.95 

4 

5% 

L2H-AMP-ST 

18.76 

17.85 

19.05 

20.00 

17.65 

16.38 

11.93 

HR-AMP-TV-2D-B 

20.88 

19.66 

20.67 

22.85 

19.60 

18.36 

16.13 



L2H-AMP-TV-2D-B 

20.47 

19.39 

20.56 

22.25 

19.10 

18.23 

16.11 


TABLE VII 

PSNRS (DB) of 128 X 128 IMAGE RECONSTRUCTIONS WITH HR-AMP AND L2H-AMP. THE TRANSFORM DOMAIN IN 

AMP-ST IS DCT. 


d=2 for LR-AMP-ST 

d=2 for LR-AMP-TV-2D 

6i% 

HR-AMP-ST 

LR-AMP-ST 

di% 

HR-AMP-TV-2D 

LR-AMP-TV-2D-R 

LR-AMP-TV-2D-B 

5 

10.8969 

3.7318 

5 

9.9753 

2.4964 

2.3032 

10 

11.5907 

3.9249 

10 

6.9049 

2.2856 

2.0812 

20 

12.5869 

4.1898 

20 

5.7327 

2.6401 

2.4869 

d=4 for LR-AMP-ST 

d=4 for LR-AMP-TV-2D 

Si% 

HR-AMP-ST 

LR-AMP-ST 

6i% 

HR-AMP-TV-2D 

LR-AMP-TV-2D-R 

LR-AMP-TV-2D-B 

3 

0.2489 

0.0075 

3 

14.4794 

0.8791 

0.8486 

4 

0.3205 

0.0080 

4 

11.8068 

0.8831 

0.8594 

5 

0.3937 

0.0107 

5 

9.9753 

0.9104 

0.9005 


TABLE VIII 

CPU RUNNING TIME IN SECONDS OE DIFFERENT METHODS FOR THE 128 X 128 BARBARA IMAGE. 


5 ) Imaging in the presence of measurement noise: Table VI shows the performance of MR-AMP in 
different domains when various amounts of measurement noises are added. The proposed LR-AMP still 
outperforms the HR-AMP and H2L-AMP in almost all cases. 

6) LR approximation: Another important problem in MR-CS is how to use a recovered LR image 
by LR-AMP to help the reconstruction of a higher-resolution image. As an initial attempt, we show in 


Table VII some results by simply upsampling the recovered LR image with the upsampling matrix to get 
a HR image, named L2H-AMP. As shown by the table, even this simple method can sometimes provide 
better HR images than HR-AMP. For example, L2H-AMP-ST can outperform HR-MP-ST in almost all 
cases. However, HR-AMP-TV-2D-B outperforms L2H-AMP-TV-2D-B when d = 4 and <5i = 0.05, which 
implies that L2H-AMP is far from optimal. The reason is that high frequency information can be captured 
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in CS measurements y, but L2H-AMP is based on LR-AMP. It thus treats the high frequency information 
as approximation errors, and the upsampling matrix cannot estimate such information from LR image. 

7) Computational complexity: The computational complexities of various methods are reported in 
Table VIII[ which shows that when d = 2, the proposed LR-AMP is about 2 times faster than the 
HR-AMP (the H2L-AMP is even slower than HR-AMP due to the additional downsampling), and the 
spatial-domain method is faster than the transform-domain one. However, when d = 4 (the size of the 
LR image is 1/16 of the HR one), the thresholding rule in soft-thresholding denoiser is changed from 


the time-consuming optimal SURE method in |30| for d = 2 to the fast suboptimal max-min method 


in | [^ . Thus, the LR-AMP-ST is about 36 times faster than HR-AMP-ST, the latter is about 25 times 
faster than the HR-AMP-TV, and LR-AMP-ST is about 100 times faster than LR-AMP-TV. Moreover, 
LR-AMP-TV is about 13 times faster than HR-AMP-TV. This gives some guidelines on how to choose 
the appropriate method according to the value of d when the complexity is a primary concern. 


VH. Conclusion and future work 

In this paper, we systematically study the multi-resolution compressed sensing reconstruction problem, 
which can stably recover a low-resolution signal when the sampling rate is too low to recover the 
full resolution signal. We develop an AMP-based solution and study its theoretical performance. We 
also develop the appropriate up-/down-sampling operators in both transform and spatial domain. The 
performance of the proposed scheme is demonstrated via simulation results. 

The proposed scheme can be further improved or applied to other applications. For example, in pOl , 
(H), the authors introduce various latest image denoising algorithms into AMP. Better performance can 
be achieved if proper up-/down-sampling matrices can be designed for these denoisers. Another topic is 
to make full use of the LR-AMP to reconstruct better HR image, i.e., to improve the performance of the 


L2H-AMP in Sec. VI-D6 Moreover, the proposed MR-AMP framework can also be applied to videos 


and multi-view images and videos |20|. 


Appendix 


In this appendix, we prove that the condition in Corollary 4.3 i.e., M(ei|r/) is a concave function of 
El, holds for the piecewise constant family in Eq. Q. 
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We start by defining a special family of distributions for simple sparse signals: 

: IE„„^ {||s[2 : ni]||o} ^ mei} , (30) 

where s[2 : ni] refers to the subvector of a signal s from the second entry to the last entry. 

Consider a signal x in the piecewise constant signal family define s as follows. 

s = [x[l],x[2] - x[l],.. .,x[ni] - x[ni - Ijf. 

If is clear that s ~ Vm where Vm G ^nfei- Therefore a bijection relationship holds between 
and since every signal generated from a distribution in paired with exactly one signal 

from every signal from paired with exactly one signal from ^ result, the 

proof in Q for the concavity of M{£i\ri) for block-sparse signals is applicable to the piecewise constant 
family. However, the proof in Q (at the end of Page 3406) was very brief. Therefore we include the 
following details for completeness. 

The goal of the concavity proof is to show that 


M{q£i -h (1 - q)e2\r]) > qM{£i\ri) -f (1 - q)M{£2\r]). 


(31) 


First, from Eq. (|^ and ([^, if a distribution vi G Tni,qei+{i-q)£2^ = ^^2 + (1 — 

where V 2 G -T'ni,ei and G because any measure in can be written as a convex 

combination of measures in and measures in Next, note that M(ei|7?) in Eq. (12i is 

obtained by tuning the denoising parameters to minimize the MSE of the least favorable distribution in 
the family. Eq. ( |^ can be proved by combining the two facts, because each term in the right hand side 
can be tuned independently to minimize its own least favorable MSE, whereas there is only one set of 
tuning parameters in the left hand side, leading to larger minimax MSE. 
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